/*
  File:    Oxygen.cpp
  Purpose: Provides routines to solve oxygen perfusion, 
           uptake, and decay. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdlib>
#include <cmath>
#include <iostream>

using namespace std;

#include "../Matrix.h"
#include "../ImprovedMath.h"
#include "DiffusionSolverSteadyState.h"
#include "Polynomial.h"

#include "NecrosisLevelSets.h"
#include "../LogfileEntry.h"

#include "../AngiogenesisInterface.h"
#include "Oxygen.h"

#include "../ParseText.h"

Polynomial o2uptakePolynomial;
Polynomial External_o2uptakePolynomial;

double External_o2polynomial_ECMtaxisLeftCutoff;
double external_o2_uptake_polynomial_ecm_left_cutoff;
double external_o2_uptake_polynomial_ecm_right_cutoff;


Polynomial o2_bulk_delivery_polynomial;
double o2_bulk_delivery_polynomial_ecm_left_cutoff;
double o2_bulk_delivery_polynomial_ecm_right_cutoff;


bool Disable_bulk_o2_source_degradation;
bool Disable_exterior_o2_uptake_reduction_with_ECM;

double Bulk_o2_source_degradation_rate;

bool CreateO2uptakeFunction( void )
{
 extern double NecrosisOxygenLevel;
 extern double QuiescenceOxygenLevel;
 extern double NecroticOxygenUptakeRate;
 extern double QuiescentOxygenUptakeRate;
 extern double TumorOxygenUptakeRate;
 
 extern char* ParameterFile; 

 Disable_bulk_o2_source_degradation = (bool) RetrieveDoubleParameter( ParameterFile, "Disable_bulk_o2_source_degradation" );
 Disable_exterior_o2_uptake_reduction_with_ECM = (bool) RetrieveDoubleParameter( ParameterFile, "Disable_exterior_o2_uptake_reduction_with_ECM" );

 Bulk_o2_source_degradation_rate = RetrieveDoubleParameter( ParameterFile , "Bulk_o2_source_degradation_rate" );

 Matrix Coeff = zeros(5,5);
 Matrix RHS = zeros(1,5);
 
 Matrix Soln = zeros(1,5);
 
 double middle_level = 0.5*( NecrosisOxygenLevel + QuiescenceOxygenLevel );
 
 *Coeff(0,0) = 1;
 *Coeff(0,1) = NecrosisOxygenLevel;
 *Coeff(0,2) = IntPow(NecrosisOxygenLevel,2);
 *Coeff(0,3) = IntPow(NecrosisOxygenLevel,3);
 *Coeff(0,4) = IntPow(NecrosisOxygenLevel,4);
 *RHS(0) = NecroticOxygenUptakeRate;
 
 *Coeff(1,0) = 0;
 *Coeff(1,1) = 1;
 *Coeff(1,2) = 2*NecrosisOxygenLevel;
 *Coeff(1,3) = 3*IntPow(NecrosisOxygenLevel,2);
 *Coeff(1,4) = 4*IntPow(NecrosisOxygenLevel,3);
 *RHS(1) = 0;
 
 *Coeff(2,0) = 1;
 *Coeff(2,1) = middle_level;
 *Coeff(2,2) = IntPow(middle_level,2);
 *Coeff(2,3) = IntPow(middle_level,3);
 *Coeff(2,4) = IntPow(middle_level,4);
 *RHS(2) = QuiescentOxygenUptakeRate;
 
 *Coeff(3,0) = 1;
 *Coeff(3,1) = QuiescenceOxygenLevel;
 *Coeff(3,2) = IntPow(QuiescenceOxygenLevel,2);
 *Coeff(3,3) = IntPow(QuiescenceOxygenLevel,3);
 *Coeff(3,4) = IntPow(QuiescenceOxygenLevel,4);
 *RHS(3) = TumorOxygenUptakeRate;
 
 *Coeff(4,0) = 0;
 *Coeff(4,1) = 1;
 *Coeff(4,2) = 2*QuiescenceOxygenLevel;
 *Coeff(4,3) = 3*IntPow(QuiescenceOxygenLevel,2);
 *Coeff(4,4) = 4*IntPow(QuiescenceOxygenLevel,3);
 *RHS(4) = 0;
 
 Soln = GaussianEliminationPartialPivoting( Coeff, RHS );
 
 o2uptakePolynomial.SetDegree(4);
 for( int i=0; i <= o2uptakePolynomial.GetDegree() ; i++ )
 { o2uptakePolynomial.SetCoefficient( i, *Soln(i) ); }
 
 
 /*
   Now do the same thing for the external one; This one satisfies:
   
   External_o2uptakePolynomial(original_ecm):
   
   External_o2uptakePolynomial( external_o2_uptake_polynomial_ecm_left_cutoff ) = 0.0;
   External_o2uptakePolynomial'( external_o2_uptake_polynomial_ecm_left_cutoff ) = 0.0;
   External_o2uptakePolynomial( external_o2_uptake_polynomial_ecm_right_cutoff ) = 1.0;
   External_o2uptakePolynomial'( external_o2_uptake_polynomial_ecm_right_cutoff ) = 0.0;
 */
 
 External_o2uptakePolynomial.SetVariable('E');

 external_o2_uptake_polynomial_ecm_left_cutoff = RetrieveDoubleParameter(ParameterFile,"external_o2_uptake_polynomial_ecm_left_cutoff");
 external_o2_uptake_polynomial_ecm_right_cutoff = RetrieveDoubleParameter(ParameterFile,"external_o2_uptake_polynomial_ecm_right_cutoff");

 External_o2uptakePolynomial.SetDegree(3);
 for( int i=0 ; i <= External_o2uptakePolynomial.GetDegree() ; i++ )
 { External_o2uptakePolynomial.SetCoefficient( i , 1.0 ); }
 
 Polynomial Prime = Derivative( External_o2uptakePolynomial );
 
 Coeff = zeros( External_o2uptakePolynomial.GetDegree()+1,External_o2uptakePolynomial.GetDegree()+1 );
 RHS = zeros(1,External_o2uptakePolynomial.GetDegree()+1 );
 
 for( int j=0; j <= External_o2uptakePolynomial.GetDegree() ; j++ )
 { *Coeff(0,j) = External_o2uptakePolynomial.GetCoefficient(j)*IntPow( external_o2_uptake_polynomial_ecm_left_cutoff , j ); }
 *RHS(0) = 0;
 *Coeff(1,0) = 0.0;
 for( int j=1; j <= External_o2uptakePolynomial.GetDegree() ; j++ )
 { *Coeff(1,j) = Prime.GetCoefficient(j-1)*IntPow( external_o2_uptake_polynomial_ecm_left_cutoff , j-1 ); }
 *RHS(1) = 0.0;
 
 for( int j=0; j <= External_o2uptakePolynomial.GetDegree() ; j++ )
 { *Coeff(2,j) = External_o2uptakePolynomial.GetCoefficient(j)*IntPow( external_o2_uptake_polynomial_ecm_right_cutoff , j ); }
 *RHS(2) = 1.0;
 for( int j=1; j <= External_o2uptakePolynomial.GetDegree() ; j++ )
 { *Coeff(3,j) = Prime.GetCoefficient(j-1)*IntPow( external_o2_uptake_polynomial_ecm_right_cutoff , j-1 ); }
 *RHS(3) = 0.0;
 
 Soln = GaussianEliminationPartialPivoting( Coeff, RHS  );
 
 for( int i=0; i <= External_o2uptakePolynomial.GetDegree() ; i++ )
 { External_o2uptakePolynomial.SetCoefficient( i , *Soln(i) ); }
 Prime = Derivative( External_o2uptakePolynomial );
 
 cout << "External o2 uptake polynomial details: " << endl 
      << "-------------------------------" << endl
	  << "p(Original_ECM) = "; 
 External_o2uptakePolynomial.display();
 cout << endl;
 
 cout << "At " << external_o2_uptake_polynomial_ecm_left_cutoff << " : p(Original_ECM) = " 
      << External_o2uptakePolynomial(external_o2_uptake_polynomial_ecm_left_cutoff) << " p'(Original_ECM) = " 
	  << Prime(external_o2_uptake_polynomial_ecm_left_cutoff) << endl;
 
 cout << "At " << external_o2_uptake_polynomial_ecm_right_cutoff << " : p(Original_ECM) = " 
      << External_o2uptakePolynomial(external_o2_uptake_polynomial_ecm_right_cutoff) << " p'(Original_ECM) = " 
	  << Prime(external_o2_uptake_polynomial_ecm_right_cutoff) << endl << endl;
	  
 /*
   Now do the same thing for the o2 bulk delivery; This one satisfies:
   
   o2_bulk_delivery_polynomial(original_ecm):
   
   o2_bulk_delivery_polynomial( o2_bulk_delivery_polynomial_ecm_left_cutoff ) = 0.0;
   o2_bulk_delivery_polynomial'( o2_bulk_delivery_polynomial_ecm_left_cutoff ) = 0.0;
   o2_bulk_delivery_polynomial( o2_bulk_delivery_polynomial_ecm_right_cutoff ) = 1.0;   
   o2_bulk_delivery_polynomial'( o2_bulk_delivery_polynomial_ecm_right_cutoff ) = 0.0;   
 */
 
 o2_bulk_delivery_polynomial.SetVariable('E');

 o2_bulk_delivery_polynomial_ecm_left_cutoff = RetrieveDoubleParameter(ParameterFile,"o2_bulk_delivery_polynomial_ecm_left_cutoff");
 o2_bulk_delivery_polynomial_ecm_right_cutoff = RetrieveDoubleParameter(ParameterFile,"o2_bulk_delivery_polynomial_ecm_right_cutoff");


 o2_bulk_delivery_polynomial.SetDegree(3);
 for( int i=0 ; i <= o2_bulk_delivery_polynomial.GetDegree() ; i++ )
 { o2_bulk_delivery_polynomial.SetCoefficient( i , 1.0 ); }
 
 Prime = Derivative( o2_bulk_delivery_polynomial );
 
 Coeff = zeros( o2_bulk_delivery_polynomial.GetDegree()+1,o2_bulk_delivery_polynomial.GetDegree()+1 );
 RHS = zeros(1,o2_bulk_delivery_polynomial.GetDegree()+1 );
 
 for( int j=0; j <= o2_bulk_delivery_polynomial.GetDegree() ; j++ )
 { *Coeff(0,j) = o2_bulk_delivery_polynomial.GetCoefficient(j)*IntPow( o2_bulk_delivery_polynomial_ecm_left_cutoff , j ); }
 *RHS(0) = 0;
 *Coeff(1,0) = 0.0;
 for( int j=1; j <= o2_bulk_delivery_polynomial.GetDegree() ; j++ )
 { *Coeff(1,j) = Prime.GetCoefficient(j-1)*IntPow( o2_bulk_delivery_polynomial_ecm_left_cutoff , j-1 ); }
 *RHS(1) = 0.0;
 
 for( int j=0; j <= o2_bulk_delivery_polynomial.GetDegree() ; j++ )
 { *Coeff(2,j) = o2_bulk_delivery_polynomial.GetCoefficient(j)*IntPow( o2_bulk_delivery_polynomial_ecm_right_cutoff , j ); }
 *RHS(2) = 1.0;
 for( int j=1; j <= o2_bulk_delivery_polynomial.GetDegree() ; j++ )
 { *Coeff(3,j) = Prime.GetCoefficient(j-1)*IntPow( o2_bulk_delivery_polynomial_ecm_right_cutoff , j-1 ); }
 *RHS(3) = 0.0; 
 
 Soln = GaussianEliminationPartialPivoting( Coeff, RHS  );
 
 for( int i=0; i <= o2_bulk_delivery_polynomial.GetDegree() ; i++ )
 { o2_bulk_delivery_polynomial.SetCoefficient( i , *Soln(i) ); }
 Prime = Derivative( o2_bulk_delivery_polynomial );
 
 cout << "o2 bulk delivery polynomial details: " << endl 
      << "-------------------------------" << endl
      << "p(Original_ECM) = "; 
 o2_bulk_delivery_polynomial.display();
 cout << endl;
 
 cout << "At " << o2_bulk_delivery_polynomial_ecm_left_cutoff << " : p(Original_ECM) = " 
      << o2_bulk_delivery_polynomial(o2_bulk_delivery_polynomial_ecm_left_cutoff) << " p'(Original_ECM) = " 
	  << Prime(o2_bulk_delivery_polynomial_ecm_left_cutoff) << endl;
 
 cout << "At " << o2_bulk_delivery_polynomial_ecm_right_cutoff << " : p(Original_ECM) = " 
      << o2_bulk_delivery_polynomial(o2_bulk_delivery_polynomial_ecm_right_cutoff) << " p'(Original_ECM) = " 
      << Prime(o2_bulk_delivery_polynomial_ecm_right_cutoff) << endl << endl;

 return true;
}

double o2uptakeFunction( double o2 )
{
 extern double NecrosisOxygenLevel;
 extern double NecroticOxygenUptakeRate;
 extern double QuiescenceOxygenLevel;
 extern double TumorOxygenUptakeRate;

 if( o2 < NecrosisOxygenLevel )
 { return NecroticOxygenUptakeRate; }
 
 if( o2 > QuiescenceOxygenLevel )
 { return TumorOxygenUptakeRate; }
 
 return o2uptakePolynomial( o2 );
}

double OxygenDiffusionConstantFunction( double u,double x, double y,int i, int j)
{
 extern Matrix OxygenDiffusivity;
 return *OxygenDiffusivity(i,j);							
}

double O2_OneFunction(double input )
{ return 1.0; }

double O2_ZeroFunction(double input )
{ return 0.0; }

double OxygenSourceTimesO2( double o2, double x, double y, int i, int j )
{
 extern double _EPS;
 
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix ECM;
 extern Matrix BloodVesselDensity;
 extern Matrix Pressure;
 extern Matrix BloodNodePressure;
 
 extern double ECMoxygenUptakeRate;
 extern double NecroticOxygenUptakeRate;
 extern double NecrosisOxygenLevel;
 extern double NecroticOxygenUptakeRate;
 
 extern double QuiescenceOxygenLevel;
 extern double QuiescentOxygenUptakeRate; 
 extern double TumorOxygenUptakeRate; 
 extern double OxygenUptakeRate;
 extern Matrix VesselCells;
 extern bool TumorBegin;
 extern bool FlowBegin;
 extern Matrix Flow;
 // Use Original_ECM instead of ECM because Original_ECM is not re-created. 
 // We should really make a new variable referring to original tissue, which is 
 // degraded simultaneously with ECM. Then, o2 bulk supply would 
 // change to bvd*original_tissue * ( bvd*original_tissue - o2 )
 // In fact, original tissue = 1 everywhere at the start? 
 // extern Matrix Original_ECM;
   extern Matrix Original_ECM;
  extern Matrix BloodVesselIndicator;
  extern double BloodTissueTransferRate;
  extern Matrix HaematocritLevel;
  //double h = *HaematocritLevel(i,j)* *Flow(i,j);
   double h = *HaematocritLevel(i,j)/0.45;
  double FlowAverage = 2.8e-16;
  
  double output = 0;
  // double f=*Flow(i,j)/FlowAverage;
  //if(f>1)
  //{
  //f=1;
  //}

if(*BloodVesselIndicator(i,j)>0.25)
{
if(!FlowBegin)
{
output += BloodTissueTransferRate;
}
else
{
double extravasation = pressure_extravasation( *Pressure(i,j) , *BloodNodePressure(i,j) );
h = h * extravasation;
output += h*BloodTissueTransferRate;
} 

}


 if( Original_ECM.TellRows() <= 1 || Original_ECM.TellCols() <= 1 )
 { Original_ECM = ECM; }
 

 if(TumorBegin)
 {
 
 if( *Phi(i,j) > _EPS )
 {
  if( *Original_ECM(i,j) <= external_o2_uptake_polynomial_ecm_left_cutoff )
  { return 0.0-output; }
  
  if( *Original_ECM(i,j) >= external_o2_uptake_polynomial_ecm_right_cutoff )
  { return -ECMoxygenUptakeRate-output; }
  
  if( Disable_exterior_o2_uptake_reduction_with_ECM )
  {
   return -ECMoxygenUptakeRate-output;
  }
  
  return -ECMoxygenUptakeRate*External_o2uptakePolynomial( *Original_ECM(i,j) )-output;
 }


 if( *PhiNecrotic(i,j) <= _EPS  )
 { return -NecroticOxygenUptakeRate-output; } 
 if( o2 <= NecrosisOxygenLevel )
 { return -NecroticOxygenUptakeRate-output; } 
 
 return -o2uptakeFunction( o2 )-output;

}
else
{
return -output -OxygenUptakeRate;
}



 
}

double OxygenWTSourceTimesO2( double o2, double x, double y, int i, int j )
{
	extern double _EPS;
	
	extern Matrix Phi;
	extern Matrix PhiNecrotic;
	extern Matrix ECM;
	extern Matrix BloodVesselDensity;
	extern Matrix Pressure;
	extern Matrix BloodNodePressure;
	extern Matrix InterPressure;
	extern Matrix HomPressure;
	extern double VFWater;
	extern double VFCell;
	extern double ChMobility;
	extern double ChLength;
	extern double MitosisRate;
	extern double ECMoxygenUptakeRate;
	extern double NecroticOxygenUptakeRate;
	extern double NecrosisOxygenLevel;
	extern double NecroticOxygenUptakeRate;
	extern double Oxygen_pi;
	extern double QuiescenceOxygenLevel;
	extern double QuiescentOxygenUptakeRate; 
	extern double TumorOxygenUptakeRate; 
	extern double OxygenUptakeRate;
	extern Matrix VesselCells;
	extern bool TumorBegin;
	extern bool FlowBegin;
	extern Matrix Flow;
	// Use Original_ECM instead of ECM because Original_ECM is not re-created. 
	// We should really make a new variable referring to original tissue, which is 
	// degraded simultaneously with ECM. Then, o2 bulk supply would 
	// change to bvd*original_tissue * ( bvd*original_tissue - o2 )
	// In fact, original tissue = 1 everywhere at the start? 
	// extern Matrix Original_ECM;
	extern Matrix Original_ECM;
	extern Matrix BloodVesselIndicator;
	extern double BloodTissueTransferRate;
	extern Matrix HaematocritLevel;
	extern Matrix HomPressure;
	//double h = *HaematocritLevel(i,j)* *Flow(i,j);
	double h = *HaematocritLevel(i,j)/0.45;
	double FlowAverage = 2.8e-16;
	/*
	if((*InterPressure(i,j)*VFWater + (MitosisRate*ChLength*ChLength* *Pressure(i,j)/ChMobility)*VFCell)>0.1*HomPressure)
	{
		h = 0;	
	}
	 */
	double output = 0;
	// double f=*Flow(i,j)/FlowAverage;
	//if(f>1)
	//{
	//f=1;
	//}
	
	if(*BloodVesselIndicator(i,j)>0.25)
	{
		if(!FlowBegin)
		{
			output += BloodTissueTransferRate;
		}
		else
		{
			
			//double extravasation = 1 - *InterPressure(i,j)/HomPressure;
			double extravasation = 1 - Oxygen_pi* *InterPressure(i,j)/ *HomPressure(i,j);
			//double extravasation = 1 - *InterPressure(i,j)*(VFWater/HomPressure) - (MitosisRate*ChLength*ChLength* *Pressure(i,j)/ChMobility)*(VFCell/HomPressure);
			if (extravasation<0)
			{
				extravasation = 0;
			}
			else if(extravasation>1)
			{
				extravasation = 1;
			}
			
			h = h * extravasation;
			output += h*BloodTissueTransferRate;
		} 
		
	}
	
	
	if( Original_ECM.TellRows() <= 1 || Original_ECM.TellCols() <= 1 )
	{ Original_ECM = ECM; }
	
	
	if(TumorBegin)
	{
		
		if( *Phi(i,j) > _EPS )
		{
			if( *Original_ECM(i,j) <= external_o2_uptake_polynomial_ecm_left_cutoff )
			{ return 0.0-output; }
			
			if( *Original_ECM(i,j) >= external_o2_uptake_polynomial_ecm_right_cutoff )
			{ return -ECMoxygenUptakeRate-output; }
			
			if( Disable_exterior_o2_uptake_reduction_with_ECM )
			{
				return -ECMoxygenUptakeRate-output;
			}
			
			return -ECMoxygenUptakeRate*External_o2uptakePolynomial( *Original_ECM(i,j) )-output;
		}
		
		
		if( *PhiNecrotic(i,j) <= _EPS  )
		{ return -NecroticOxygenUptakeRate-output; } 
		if( o2 <= NecrosisOxygenLevel )
		{ return -NecroticOxygenUptakeRate-output; } 
		
		return -o2uptakeFunction( o2 )-output;
		
	}
	else
	{
		return -output -OxygenUptakeRate;
	}
	
	
	
	
}

 double OxygenWTSourceRemainder( double o2, double x, double y, int i, int j )
{ 
	extern Matrix PhiNecrotic;
	extern Matrix PhiQuiescent;
	extern Matrix Phi;
	extern double _EPS;
	extern double dx;
	extern double dy;
	extern double BloodTissueTransferRate;
	extern Matrix BloodVesselIndicator;
	extern bool TumorBegin;
	extern bool FlowBegin;
	extern Matrix HaematocritLevel;
	extern Matrix Flow;
	double HaematocritScale = 0.5; 
	double HaematocritToOxygenScale = 10;
	extern Matrix VesselCells;
	extern Matrix InterPressure;
	extern Matrix Pressure;
	extern Matrix HomPressure;
	extern double Oxygen_pi;
	extern double VFWater;
	extern double VFCell;
	extern double ChMobility;
	extern double ChLength;
	extern double MitosisRate;

	extern int Xnodes,Ynodes;
	double FlowRef=1.9096e-11; //m3/s (Flow in PV)
	double FlowAverage = 2.8e-16;
	extern Matrix BloodNodePressure;
	
	extern Matrix BloodVesselPressureX;
	extern Matrix BloodVesselPressureY;
	
	if( HaematocritLevel.TellRows() != Xnodes ||
	   HaematocritLevel.TellCols() != Ynodes )
	{ HaematocritLevel = zeros(Xnodes,Ynodes); }

	
	double h = *HaematocritLevel(i,j)/0.45;
	/*
	if((*InterPressure(i,j)*VFWater + (MitosisRate*ChLength*ChLength* *Pressure(i,j)/ChMobility)*VFCell)>0.1*HomPressure)
	{
		h = 0;	
	}
	 */
	
	
	//double extravasation = 1 - *InterPressure(i,j)/HomPressure;
	double extravasation = 1 - Oxygen_pi* *InterPressure(i,j)/ *HomPressure(i,j);
	//double extravasation = 1 - *InterPressure(i,j)*(VFWater/HomPressure) - (MitosisRate*ChLength*ChLength* *Pressure(i,j)/ChMobility)*(VFCell/HomPressure);
	if (extravasation<0)
	{
		extravasation = 0;
	}
	else if(extravasation>1)
	{
		extravasation = 1;
	}
	

	h = h * extravasation;
	
	
	
	if(*BloodVesselIndicator(i,j)>0.25)
	{
		if(!FlowBegin)
		{
			return BloodTissueTransferRate;
		}
		else
		{
			//return f*BloodTissueTransferRate;
			return h*BloodTissueTransferRate;
		} 
		
		
	}
	else
	{
		return 0;
	}
}

double OxygenSourceRemainder( double o2, double x, double y, int i, int j )
{
 extern Matrix Phi;
 extern Matrix Structure;

 extern double bulk_blood_tissue_transfer_rate;
 extern double bulk_blood_transfer_cutoff_distance;
 extern double MinimumMobility;
 extern Matrix BloodVesselDensity;
 extern Matrix ECM;

 double output = 0.0;
 
 extern double _EPS;
 
 
 /*
 // 6-13-2007
    Disable_bulk_o2_source_degradation = (bool) RetrieveDoubleParameter( ParameterFile, "Disable_bulk_o2_source_degradation" );
 Disable_exterior_o2_uptake_reduction_with_ECM = (bool) RetrieveDoubleParameter( ParameterFile, "Disable_exterior_o2_uptake_reduction_with_ECM" );
*/
 
 
/*
 // 6-25-2007
 double b = *BloodVesselDensity(i,j);
 output = bulk_blood_tissue_transfer_rate*b*(b-o2); 
 if( output < 0.0 || *Phi(i,j) <= _EPS || *Structure(i,j) > 0.1 )
 { output = 0.0; }
*/ 
 
 double b = *BloodVesselDensity(i,j);
// extern Matrix Original_ECM;
 //double original_ecm = *Original_ECM(i,j);
 
 /*
 // disabled on 6-27-2007
 double poly = o2_bulk_delivery_polynomial( original_ecm );
 if( original_ecm < o2_bulk_delivery_polynomial_ecm_left_cutoff )
 { poly = 0.0; }
 if( original_ecm > o2_bulk_delivery_polynomial_ecm_right_cutoff )
 { poly = 1.0; }
 */
 double poly = 1.0;
 /*
 output = bulk_blood_tissue_transfer_rate*poly*b*(1.0-o2); 
 if( output < 0.0 || *Phi(i,j) <= _EPS || *Structure(i,j) > 0.1 )
 { output = 0.0; }
 */

/* 
 double ecm = *ECM(i,j); 

 if( *Phi(i,j) > _EPS && 
     ecm >= 0.2 && 
     *Structure(i,j) < 0.1 && 
     *BloodVesselDensity(i,j) > 0.01 )
 {
//  double b = *BloodVesselDensity(i,j);
 
  output = bulk_blood_tissue_transfer_rate*(1.0-o2)*b; 
 }
*/ 
 
 return output + TransferOxygenFromBloodVessels( o2,0.0,i,j );
}

double EmptyFunction(double input)
{ return 0.0; }

bool UpdateOxygenDiffusivity( void )
{
 extern double _EPS;
 
 extern int Xnodes;
 extern int Ynodes;
 extern Matrix Structure;
 
 extern Matrix OxygenDiffusivity;
 
 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0 ; j < Ynodes ; j++ )
  {
   // UNCOUPLE if( *Structure(i,j) > 0.1 )
   {
    *OxygenDiffusivity(i,j) = 1;
   }
  }
 }

 return true; 
}

bool UpdateOxygen( void )
{
 extern int Xnodes;
 extern int Ynodes;
 extern double dx;
 
 extern Matrix X;
 extern Matrix Y;
 
 extern Matrix Oxygen;
 extern Matrix OxygenDiffusivity;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 extern Matrix BloodVesselIndicator;

 extern char* StringBuffer;
 extern char* LogFile;

 extern bool TopBloodVessel;
 extern bool IsFirstOxygenUpdate;
 
 // update the diffusivity

 
// UpdateOxygenDiffusivity();


 // solve the oxygen perfusion equation 
 // to steady state;
 
 int MaxIterations;
 double Tolerance;

 if( Oxygen.TellRows() != Xnodes || Oxygen.TellCols() != Ynodes )
 { 
  Oxygen = zeros(Xnodes,Ynodes); 
  StartStringBuffer( "Warning: Oxygen previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( OxygenDiffusivity.TellRows() != Xnodes || OxygenDiffusivity.TellCols() != Ynodes )
 { 
  OxygenDiffusivity = ones(Xnodes,Ynodes);  
  StartStringBuffer( "Warning: Oxygen diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 
 /* UNCOUPLE
 if( PhiNecrotic.TellRows() != Xnodes || PhiNecrotic.TellCols() != Ynodes )
 {
  PhiNecrotic = ones(Xnodes,Ynodes);
  StartStringBuffer( "Warning: PhiN diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( PhiQuiescent.TellRows() != Xnodes || PhiQuiescent.TellCols() != Ynodes )
 {
  PhiQuiescent = ones(Xnodes,Ynodes); 
  StartStringBuffer( "Warning: PhiQ diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 */
 if( BloodVesselIndicator.TellRows() != Xnodes || BloodVesselIndicator.TellCols() != Ynodes )
 {
  BloodVesselIndicator = zeros(Xnodes,Ynodes);
  StartStringBuffer( "Warning: Blood Vessel Indicator previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }

 // test line
// UNCOUPLE PhiQuiescent = ones(Xnodes,Ynodes);
 // end test line
  
// char TopO2conditionType = 'D';
 //double (*TopO2BC)(double) = O2_OneFunction;
/* UNCOUPLE
 if( !TopBloodVessel )
 {
  TopO2conditionType = 'N'; 
  TopO2BC = O2_ZeroFunction;
 }
 
 extern bool BottomBloodVessel;
 char BottomO2conditionType = 'D';
 double (*BottomO2BC)(double) = O2_OneFunction;
 if( !BottomBloodVessel )
 {
  BottomO2conditionType = 'N'; 
  BottomO2BC = O2_ZeroFunction;
 }
   
*/
    MaxIterations = 20000;  // 20000; 
  Tolerance = IntPow(dx,3); // IntPow(dx,3);

  Oxygen = NonlinearDiffusionSteadyState(X,Y,Oxygen, 
    O2_ZeroFunction, 'N', O2_ZeroFunction, 'N', 
    O2_ZeroFunction, 'N', O2_ZeroFunction, 'N',  
    OxygenDiffusionConstantFunction,
    OxygenSourceTimesO2,
    OxygenSourceRemainder,
    MaxIterations, Tolerance );
  IsFirstOxygenUpdate = false; 
 //   cout<<"1"<<endl;

  MaxIterations = 20000;  // 20000; 
  Tolerance = IntPow(dx,4); // IntPow(dx,4);

  Oxygen = NonlinearDiffusionSteadyState(X,Y,Oxygen, 
    O2_ZeroFunction, 'N', O2_ZeroFunction, 'N', 
    O2_ZeroFunction, 'N', O2_ZeroFunction, 'N',  
    OxygenDiffusionConstantFunction,
    OxygenSourceTimesO2,
    OxygenSourceRemainder,
    MaxIterations, Tolerance );
  IsFirstOxygenUpdate = false; 
//  cout<<"1"<<endl;

/*  Oxygen = NonlinearDiffusionSteadyStateNonAdaptive(X,Y,Oxygen, 
    O2_ZeroFunction, 'N', O2_ZeroFunction, 'N', 
    BottomO2BC, BottomO2conditionType, TopO2BC, TopO2conditionType,  
    OxygenDiffusionConstantFunction,
    OxygenSourceTimesO2,
    OxygenSourceRemainder,
    MaxIterations, Tolerance );
  IsFirstOxygenUpdate = false; 
  */
//  cout<<"1"<<endl; 
  
// UpdateNecrosisLevelset();

 //UpdateQuiescentLevelset();	
 	   

 return true;
}

bool QuicklyUpdateOxygen( void )
{
 extern int Xnodes;
 extern int Ynodes;
 extern double dx;
 
 extern Matrix X;
 extern Matrix Y;
 
 extern Matrix Oxygen;
 extern Matrix OxygenDiffusivity;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 extern Matrix BloodVesselIndicator;

 extern char* StringBuffer;
 extern char* LogFile;

 extern bool TopBloodVessel;
 extern bool IsFirstOxygenUpdate;
 
 // update the diffusivity
 
 UpdateOxygenDiffusivity();

 // solve the oxygen perfusion equation 
 // to steady state;
 
 int MaxIterations;
 double Tolerance;

 if( Oxygen.TellRows() != Xnodes || Oxygen.TellCols() != Ynodes )
 { 
  Oxygen = zeros(Xnodes,Ynodes); 
  StartStringBuffer( "Warning: Oxygen previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( OxygenDiffusivity.TellRows() != Xnodes || OxygenDiffusivity.TellCols() != Ynodes )
 { 
  OxygenDiffusivity = ones(Xnodes,Ynodes);  
  StartStringBuffer( "Warning: Oxygen diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( Oxygen.TellRows() != Xnodes || Oxygen.TellCols() != Ynodes )
 {
  Oxygen = zeros(Xnodes,Ynodes);
  StartStringBuffer( "Warning: Oxygen diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( PhiNecrotic.TellRows() != Xnodes || PhiNecrotic.TellCols() != Ynodes )
 {
  PhiNecrotic = ones(Xnodes,Ynodes);
  StartStringBuffer( "Warning: PhiN diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( PhiQuiescent.TellRows() != Xnodes || PhiQuiescent.TellCols() != Ynodes )
 {
  PhiQuiescent = ones(Xnodes,Ynodes); 
  StartStringBuffer( "Warning: PhiQ diffusivity previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 if( BloodVesselIndicator.TellRows() != Xnodes || BloodVesselIndicator.TellCols() != Ynodes )
 {
  BloodVesselIndicator = zeros(Xnodes,Ynodes);
  StartStringBuffer( "Warning: Blood Vessel Indicator previously uninitialized. Building ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 }
 
 // test line
 PhiQuiescent = ones(Xnodes,Ynodes);
 // end test line
 
 char TopO2conditionType = 'D';
 double (*TopO2BC)(double) = O2_OneFunction;
 if( !TopBloodVessel )
 {
  TopO2conditionType = 'N'; 
  TopO2BC = O2_ZeroFunction;
 }
  
  MaxIterations = 2000;  // 20000; 
  Tolerance = IntPow(dx,4); // IntPow(dx,3);
    
  Oxygen = NonlinearDiffusionSteadyState(X,Y,Oxygen, 
    O2_ZeroFunction, 'N', O2_ZeroFunction, 'N', 
    O2_ZeroFunction, 'N', TopO2BC, TopO2conditionType,  
    OxygenDiffusionConstantFunction,
    OxygenSourceTimesO2,
    OxygenSourceRemainder,
    MaxIterations, Tolerance );

 return true;
}
